Force Control of a Haptic Flexible-Link Antenna Based on a Lumped-Mass Model

Haptic organs are common in nature and help animals to navigate environments where vision is not possible. Insects often use slender, lightweight, and flexible links as sensing antennae. These antennae have a muscle-endowed base that changes their orientation and an organ that senses the applied force and moment, enabling active sensing. Sensing antennae detect obstacles through contact during motion and even recognize objects. They can also push obstacles. In all these tasks, force control of the antenna is crucial. The objective of our research is to develop a haptic robotic system based on a sensing antenna, consisting of a very lightweight and slender flexible rod. In this context, the work presented here focuses on the force control of this device. To achieve this, (a) we develop a dynamic model of the antenna that moves under gravity and maintains point contact with an object, based on lumped-mass discretization of the rod; (b) we prove the robust stability property of the closed-loop system using the Routh stability criterion; and (c) based on this property, we design a robust force control system that performs efficiently regardless of the contact point with the object. We built a mechanical device replicating this sensing organ. It is a flexible link connected at one end to a 3D force–torque sensor, which is attached to a mechanical structure with two DC motors, providing azimuthal and elevation movements to the antenna. Our experiments in contact situations demonstrate the effectiveness of our control method.


Introduction
Haptics is the technology of touch.In recent years, there has been increasing interest in developing integrated sensory systems, particularly those involving various tactile sensors [1].Multiple applications require integrated systems, such as machine assemblies for precise positioning, impact protection, navigation, etc. Tactile/touch sensing is essential for developing human-machine interfaces and electronic skins in areas such as automation, security, and medical care [2].Tactile sensors were first explored in the early 1990s, for example, in the work of Russell [3].Since then, natural tactile sensors, including whiskers and antennae, have been investigated, e.g., [4].Several attempts have been made to build biomimetic active sensory applications, also known as vibrational systems.Mammal and insect sensing (see Figure 1) has inspired multiple engineering applications, such as the whisker-based texture discrimination presented in [5].The less frequent use of tactile sensing may be partly attributed to its complex and distributed nature.Issues such as sensor placement, robustness, and wiring complexity, among others, make its effective utilization challenging.
Previous works have provided evidence that artificial vibrissal systems can compute estimates of distance and shape [6][7][8] and can distinguish between textures with different spatial frequencies [9,10].These results demonstrate the potential for vibrissal sensors as effective devices for tactile object recognition.Sensors used in these systems include electret microphones [11], resistive arrays [12], strain gauges [7], piezoelectric sensors [9], and magnetic Hall-effect sensors [9,13].Each of these technologies has its advantages and disadvantages.In the last two decades, a robust and compact sensor device called the "sensing antenna" has been proposed, efficiently addressing some of the aforementioned problems.This active sensor consists of a flexible beam moved by servo-controlled motors and a load cell placed between the beam and the motors.An example of this device is shown in Figure 2. The sensing antenna replicates the touch sensors found in some animals and employs an active sensing strategy.The servomotor system moves the beam back and forth until it hits an object.At this instant, information from the motor angles combined with force and torque measurements allows us to calculate the positions of the hit points, which represent valuable information about the object surface.Using this device, a 3D map of the object's surface, which enables recognition, can be obtained [13].Recognition is carried out using techniques that combines information on partial views to gather comprehensive information about the object, e.g., [14,15].
Two strategies can be applied to obtain such 3D maps.The first strategy involves continuously moving the beam back and forth to hit the object at different points, determining their 3D coordinates and producing a map of the object surface.This strategy is used by some insects that employ their antennae for this purpose (e.g., [16,17]).The second strategy involves sliding the beam across the object while exerting a controlled force on the surface of the object to maintain contact, collecting 3D coordinates of points on the object's surface during this movement.This strategy is utilized by some mammals with whiskers as sensors (e.g., [18,19]).Both strategies can be implemented using the aforementioned sensing antennae and require precise control of the force exerted by the antenna on the object (e.g., [20]).Additionally, if the object or obstacle needs to be removed, force control is also necessary.
However, multiple constraints limit the performance of these devices, such as their flexible beam length, light weight, and flexibility.These characteristics make the dynamic behavior of these antennae exhibit an infinite number of vibration modes, resulting in dynamic models of infinite dimension [21].This complexity makes it very difficult to accurately control the position of these devices or achieve precise force control.If the control system in charge of moving the motors does not consider these dynamics, i.e., the beam elasticity, residual vibrations appear that prevent the accurate and fast achievement of the desired pushing force on the object being searched or moved.Moreover, permanent collisions with the object could occur, where the antenna continuously moves back and forth as it collides with the object.These phenomena would cause delays in the recognition process and diminish the quality of object surface estimates, and therefore reduce the efficiency of the device's functioning.
It is well known that the amplitudes of the vibration modes of flexible beams diminish as the frequencies of the modes increase.This allows us to truncate their infinite-order models to finite-order models (e.g., [22]) that usually include as many as three or four vibration modes, which yield accurate approximations of the antennae dynamics.This truncation can be applied to the dynamic model obtained using an assumed-modes modeling approach (e.g., [23]), or directly to the mechanism by assuming lumped masses (e.g., [24]).
Assumed modes (e.g., [23]) have already been used to model the contact dynamics of a beam (or flexible antenna) [25] in the horizontal plane, but never in the vertical plane, under the effect of gravity.Moreover, these models are complex and pose difficulties in designing robust control systems for the contact force.Lumped mass models have been developed for the free rotation movements of flexible beams in the horizontal plane, in attitude movements where gravity affects the dynamics [26], or in a two-degrees-offreedom single flexible-link antenna [27].Robust force control at the tip of the beam has been addressed using a lumped-mass model with a single lumped mass at the tip in [28] in the case of having a horizontal degree of freedom.However, a lumped-mass model has never been used to model the dynamics of the contact situation in which a flexible beam pushes an object at an intermediate point along its link.This requires the use of several lumped masses.
The control of the force at the tip of a single flexible link that rotates on an horizontal plane through one of its ends was studied in [29], assuming a distributed mass link.The experiments showed that direct force feedback from a sensor placed at the tip could not ensure closed-loop stability.Stable tip contact control of a distributed mass link, where a switching transition occurred between the unconstrained and constrained environments, was achieved by [30] using a PD controller that provides feedback from hub measurements.This yielded a control system robust to the mechanical impedance of the contacted object, but could not achieve force control.To increase the stability of the tip force control, some works have redefined the force output to be fed back, e.g., [31,32].In [33], the tip-contactforce control of a constrained single-link flexible arm was performed, overcoming the non-minimum phase nature of the system by defining a new input and generating a virtual contact-force output through a parallel compensator.It was proven that the transfer function from the new input to the virtual contact-force output was minimum-phase and stable.Ref. [34] also addresses tip-contact-force control of a one-link flexible arm interacting with a rigid environment.To achieve contact-force control, a boundary controller was proposed based on an infinite-dimensional dynamic model.The contact-force control and vibration suppression problem for a constrained one-link flexible manipulator with an unknown control direction and a time-varying actuator fault was studied in [35].Finally, we again mention [28], where fractional-order control was implemented in a massless link with a tip payload, damping rebounds and ensuring robust stability to the mechanical impedance of the contacted object.
All these works focus on the force control of a flexible link interacting with the environment, considering that contact is made at the tip.To the best of our knowledge, the force control at an intermediate point of a flexible beam in an elevation rotation movement has never been addressed, either using an assumed-modes model or a lumped-mass model.The objectives of the present research are as follows: (1) to establish a model of the dynamics of a flexible beam contacting an object at one of its intermediate points in a rotational elevation movement based on multiple lumped masses, and (2) based on that model, to define a control system that exerts a precise pushing force on an object.Both objectives represent the contributions of this paper and have never been previously addressed.
This paper is organized is as follows.After this Introduction, Section 2 presents our experimental setup of a flexible-link antenna.Section 3 develops the dynamic model based on lumped masses.Section 4 fits this model to the lowest-frequency mode obtained from an assumed-mass model.Section 5 obtains the transfer functions of our prototype, and Section 6 derives a robust control system based on these functions.Section 7 presents our experimental results, and Section 8 offers our conclusions.

Experimental Setup
The experimental prototype is a two-degrees-of-freedom (2DOF) robotic system with a single flexible link, which is used as a sensing antenna in haptics applications.A detailed 3D representation is shown in Figure 2. Its design was developed by our group in previous works [36], where it has been employed as a tactile sensor to detect objects in its surroundings.
The flexible link, also referred to as antenna, is a lightweight, slender carbon-fiber rod with a circular cross-section.It is fixed at one of its ends (the base), while the other end moves freely (the tip).The antenna is attached at the base to a six-axis ATI FTD-MINI40 force-torque (F-T) sensor, which measures the Cartesian reacting forces and torques generated by the link.The signals are acquired through gauges located inside the sensor, which are multiplexed and amplified to send the information regarding forces and torques to a data acquisition card (DAQ).Holding the sensor and the antenna there is the servomotor structure, which is driven by two Harmonic Drive PMA-5A direct-current (DC) mini-servo actuator motor sets, featuring zero-backlash 1:100 reduction gears.One servomotor rotates the system with azimuthal movements (horizontal plane), while the other rotates it with elevation movements (vertical plane).These DC motors have incremental optical encoders that measure the angular position of the motors, θ m 1 and θ m 2 , corresponding to the azimuthal and elevation joints, respectively.Additionally, a stainless-steel structure holds all this equipment and fixes the system to a flat surface with three legs to ensure perfect stability.The robot is connected to a PC through data acquisition cards.The data acquisition and control algorithms were programmed using LabVIEW NXG 5.1 with a sampling time of T s = 1 ms.All work related to data analysis and representation was carried out using MATLAB 8.2.0.29 (R2013b).

Dynamic Model
This section focuses on the modeling of a flexible beam connected to a motor.This model enables us to characterize the active sensor mentioned earlier, which moves in a vertical rotation back and forth within a plane until it makes contact with an object.Thus, the effect of gravity is considered.For this study, we assume that the interaction between the structure and the environment occurs at a single point of contact along the beam.Additionally, we assume that the force applied by the object on the beam is perpendicular to it.This assumption neglects any slipping that may occur between the two bodies.Furthermore, the contacted object is assumed to be rigid.
The dynamic model of the system is divided into two parts to describe the behavior of the motor and the flexible beam.These subsystems are interconnected through the motor angle and the torque exerted on the motor by the beam, known as the coupling torque.

Beam Dynamics
The flexible beam is characterized by its length L, linear mass density ρ, and flexural rigidity EI.It is assumed that the beam is described by a massless link with n lumped masses along its length, as presented in [24].The beam exhibits small deflections, i.e, deflection lower that 10% of L, allowing us to use a linear deflection model [37].Furthermore, the internal and external friction effects of the beam are neglected.
The deflection, denoted as z(x, t), is measured relative to its undeformed position, defined by the frame (X, Z).As illustrated in Figure 3, the frame (X, Z) rotates relative to a fixed frame (X 0 , Z 0 ).This rotation is given by the angle of the motor θ m (t).Furthermore, the lumped mass m j is located at distance l µ j and angle θ µ j (t) with respect to the axis X and the frame (X 0 , Z 0 ), respectively.It should be noted that the mass m n is placed at the tip of the beam, i.e., l µ n = L.
Non-rigid contact is defined by two angles with respect to (X 0 , Z 0 ): the equilibrium angle of the surface of the contacted object, θ e , and the angle at which the link has penetrated into the object, θ c (t).If the contact is rigid, then θ c (t) = θ e .The contact position is defined by the distance l c along the X-axis.
Considering gravity, its direction opposes the Z 0 -axis.The effect of gravity on the beam is assumed to be the force computed on the undeformed beam.This is based on the principle of superposition, which can be applied when considering small deformations [38].Therefore, the deflection of a massless beam is given by and is related to the angles of the system by means of where θ(x, t) is the angle between any point of the beam and the frame (X 0 , Z 0 ).The solution of Equation ( 1) is given by a piecewise function defined as for the interval [l i−1 , l i ] with i = 1, 2, . . ., N and with l 0 = 0 and l N = L. Here, N can be either N = n or N = n + 1 depending on whether contact occurs or at which position it occurs.The distance l i is determined by the position of either one of the lumped masses l µ j or the contact point l c .
The polynomial coefficients u i,j (t) are obtained from the following conditions where ( 4) and ( 5) are the boundary conditions at the joint with the motor and the tip of the beam, respectively.Equation ( 6) represents the continuity conditions with i = 2, . . ., N.
The force F i (t) is defined by where k is the stiffness of the contacted object.If the contact is rigid, then k tends to infinity.Additionally, as previously mentioned, gravity is applied with respect to the undeformed position of the beam.The coupling torque between the beam and the motor is given by Hereinafter, we work with the nondimensional model to ensure generality and applicability to any slewing flexible beam.Defining T = ρL 4 EI , we obtain the nondimensional time τ = t/T and frequency ω = ω d T (letting ω d be the natural frequency of the flexible beam).The nondimensional spatial coordinate and deflection are χ = x/L and ζ(χ, τ) = z(x, t)/L.The forces and their positions are defined as  i (τ) = F i (t) T 2 ρL 2 and λ i = l i /L, respectively, and the nondimensional torque is Γ(τ) = Γ(t) T 2 ρL 3 .Moreover, the masses are defined as µ i = m i ρL , gravity as ĝ = g T 2 L , and the angles as θ i (τ) = θ i (t).Thus, the nondimensional form of Equation ( 1) is where, from now on, ( ˙) and ( ′ ) denote the derivatives with respect to the nondimensional time and spatial variables, respectively.The relation between the deflection and the angles ( 2) is now and the solution of the deflection (3) is (11)   where Furthermore, the conditions (4)-( 6) and the forces (7) are The coefficients υ i,j (τ) obtained by conditions ( 12)-( 14) are presented below, while their derivation is detailed in Appendix A. For the two first coefficients, it is found that υ 1,0 (τ) = 0, υ 1,1 (τ) = 0 and with i = 2, . . ., N, whereas, the others are with i = 1, . . ., N. Therefore, using the solution of the deflection in (10), the following equation is obtained Furthermore, the coupling torque of Equation ( 8) becomes Finally, the general solution defined in (20) can be employed to derive the dynamic equations of the beam in two cases: when the beam is freely vibrating and when it is in contact with an object.

Free-Vibration Model
When the link vibrates freely, the forces are only due to the concentrated masses of the model, which can be expressed as  i (τ) = µ j λ µ i θi (τ) + ĝ cos(θ m (τ)) .Therefore, for this case, the number of intervals into which the displacement ζ(χ, τ) is divided is equal to the number of masses, N = n, and the distances λ i and angles θ i (τ) are λ µ j and θ µ j (with i, j = 1, . . ., n).
Thus, Equation (20) and the coupling torque of ( 21) become These equations are expressed in a compact form as where 1 n is a vector of ones belonging to ℜ n×1 , and Finally, by manipulating the above equations, we obtain a dynamic model for the case of free vibrations:

Contact Model
Upon establishing contact, the beam oscillates around the position of the object due to the assumption of small displacements.For this reason, the following incremental angles are defined: Furthermore, the shape and dimensions of the model will depend on the relative position between the contact and the masses.Contact may occur at an intermediate position between two masses or coincide with one of them.
In the case of contact between masses, the displacement is divided into N = n + 1 intervals.The distances λ i and angles θ i (τ) (with i = 1, . . ., N) take the values of λ µ j and θ µ j (with j = 1, . . ., n) or λ c and θ c .The equations derived from Equation ( 20) are when when λ µ j−1 < λ i = λ c < λ µ j , and when The coupling torque (21) becomes The system of equations can be expressed in a more concise form as follows: where 0 n is a vector of zeros in ℜ n×1 ; 1 n+1 is a vector of ones in ℜ n+1×1 ; and M, H, Λ 1 , and Λ 2 are the same as in ( 27)-( 30) and After calculating Equations ( 38) and ( 39), we obtain the following dynamic model: where and By assuming that the contact is rigid, i.e., ∆θ c (τ) = 0, we can derive that the Equation ( 43) is the coupling torque ( 46) is and the contact force can be obtained from (44) as On the other hand, when the contact occurs in the position of one of the masses, the displacement is divided into N = n intervals and the equations derived from Equation ( 20) are when when when Here, the coupling torque is equal to (37).
Once more, the model is expressed in a concise form as follows: where the mass at which contact occurs is designated as and Ĥc2 ∈ ℜ 1×n−1 are derived by removing from ∆θ µ , Λ 1 , Λ 2 , H c1 and H c2 the element related to the mass coincident with the contact µ c .Similarly, the matrices M ∈ ℜ n−1×n−1 and Ĥ are derived by removing from M and H the columns and rows corresponding to the mass µ c .
The following dynamic model is obtained from (53) and ( 54) where and Finally, by assuming rigid contact in Equations ( 55) and (58), we obtain and Equation (56) provides us with the contact force as follows: In summary, the contact model is described by Equations ( 47)-( 49) when the contact is between the masses, and by ( 59)-(61) when it coincides with the position of a mass.
Remark.For a given number of masses n, the rigid contact model has one order fewer when the contact is at one mass (Equations ( 59)-( 61)) compared to when it is between masses (Equations ( 47)-( 49)).Consequently, when contact occurs at one of the masses, the number of vibration frequencies is reduced by one.

Motors Dynamics
The behavior of the motor is described by the following equation where J 0 is the rotational inertia of the motor and Γ m (t) is the torque produced by the actuator to move the system.The nondimensional form of Equation ( 62) is where the nondimensional inertia of the actuator is obtained by calculating , where In the case of contact, where the incremental angles of Equation ( 33) are defined, the equation of the engine is such that

Adjustment of the Rigid Contact Model
In this section, the parameters of the lumped-mass model described by the Equations ( 47)-( 49) and ( 59)-(61) are adjusted.The aim is to ensure that the first vibration frequency of the model coincides with the frequency of a flexible beam in contact with a rigid object.This frequency, as discussed in [25], depends on the point at which the contact occurs, denoted as λ c .Furthermore, the order of the model should be kept as low as possible in order to minimize the computational complexity.It is important to note that since the model is dimensionless, the results obtained here are applicable to any slewing flexible beam.
To fit the model, we start with a lower-order model where only one mass is considered, and then increase the number of masses n until a satisfactory result is achieved for the system frequency.The parameters to be adjusted are the masses µ j and their respective distances λ µ j .The conditions imposed include keeping the total mass of the link and ensuring that the length of the link is maintained so that the last mass µ n is positioned at the end of the link, i.e., λ µ n = 1.
Remark.It should be noted that the parameters of the model with one mass (n = 1) are determined by the imposed conditions.Consequently, if these conditions are to be maintained, it is not possible to modify the parameters in order to obtain a better fit.
Thus, Matlab was used to adjust the above parameters with the aim of minimizing the following mean square error: where ω1 (λ c,i ) and ω1 (λ c,i ) represent the frequencies obtained from the model in [25] and from the lower eigenvalue of the matrix R (or R), respectively.Both frequencies depend on the contact point λ c .The contact point is defined in the interval (0, 1] with a step of 0.01, resulting in a total number of points of N p = 100.The smallest number of masses that achieves a satisfactory fit of the first frequency is n = 3.The results for the models with a lower number of masses are presented in Appendix B. To achieve the fit for the three-mass model, it is considered that µ 1 varies within the interval (0, 0.99), and µ 2 within the interval (0, µ 1 ), both with a step of 0.01.For the positions of the masses, λ µ 2 and λ µ 1 vary within the intervals (0.01, 1) and (0, λ 2 ), respectively, with a step of 0.01.Moreover, µ 3 is defined by (65) and λ µ 3 = 1.The model with the minimum mean square error, i.e., MSE = 0.008, is characterized by the parameters µ 1 = 0.49, µ 2 = 0.41, µ 3 = 0.10, λ µ 1 = 0.26, and λ µ 2 = 0.70.The comparison between ω1 (λ c ) and ω1 (λ c ) for this model is presented in Figure 4. ) to ω1 (λ c ).

Flexible Beam Transfer Functions
The transfer function that relates the coupling torque to the motor angle in rigid contact is derived using the parameters from the previous section and Equations ( 47) and (48), or Equations ( 59) and (60) if the contact coincides with one of the masses.However, for the sake of clarity, the steps will be outlined for the first set of equations, but are exactly the same for the second.
As observed, the dynamic model obtained is nonlinear because gravity depends on the motor angle.Consequently, to obtain a linear model, the equations of the system are linearized around the point ∆θ m,0 = 0 (67) Thus, taking into account that cos(θ m (τ)) = cos(∆θ m (τ) + θ e ), the equilibrium point is defined by Let us define the variations from the equilibrium points as δθ m (τ) = ∆θ m (τ) − ∆θ m,0 , δθ µ (τ) = ∆θ µ (τ) − ∆θ µ,0 and δΓ coup (τ) = Γ coup (τ) − Γ coup,0 .The following linearized model is obtained by using the first-order Taylor series expansion Taking the Laplace transforms of these equations and substituting δθ µ from (70) into (71), we obtain the transfer function between the coupling torque and the angle of the motor with Truncating the transfer function to the first mode of vibration, we obtain a function of the form where the coefficients a(λ c ), b(λ c ), c(λ c ), and d(λ c ) are obtained with Matlab and are fitted by means of the functions in the Appendix C. It is important to note that the coefficient a(λ c ) corresponds to the first vibration frequency.These coefficients are obtained for the nondimensional model.Consequently, in order to obtain a valid transfer function for the dimensional model, the following transformations must be made: and the transfer function becomes Upon calculating Equation (77), a different representation of the transfer function is obtained: where K a (λ c ) is the gain of the model, β(λ c , θ e ) represents the zeros, and α(λ c ) denotes the poles, obtained from: Evaluating the values of α(λ c ) and β(λ c , θ e ) across the entire range λ c ∈ [0, 1] and This can be verified in Figure 5, which illustrates the evolution of α(λ c ) and β(λ c , θ e ) with respect to λ c .This property is crucial for demonstrating the robustness of our force controller in the subsequent section.Note that the value of β(λ c , θ e ) barely changes with respect to θ e .This is because c * (λ c ) is significantly smaller than b * (λ c ), resulting in a β(λ c , θ e ) ≈ β(λ c ).This observation can be checked in Figure A2 in Appendix C.

Control System
This section introduces a force control system for a single-link flexible robot operating under gravity.The objective is to regulate the force exerted by the flexible link on the environment, regardless of the contact point on the link.The control process is divided into three stages: (1) free motion control, where the link is servo-controlled until it makes contact with an object; (2) post-impact, where the link pushes against the object, gathering data from the force-torque sensor and using an estimator to identify where the contact point has been produced; and (3) force control, where the force exerted on the object is regulated using the information from the previous estimator.The system utilizes feedback from measurements of the motor's position and force-torque at the base of the link.
Thus, the subsystems comprising the control system can be classified into two categories: controllers and estimators.

Controllers: (a)
The motor position controller, which comprises the inner loop, regulates the dynamics of the motors between the motor angle θ m (t) and its reference θ * m (t).This design ensures robustness against Coulomb friction, viscous friction, variations in link parameters, and external forces exerted on the link, allowing the system dynamics to be treated as a linear time-invariant system.(b) The force controller, comprising the outer loop, regulates the force applied by the antenna at the contact point to a desired value F * (t).This controller operates in conjunction with the inner loop.Thus, this outer loop uses feedback measurements of force-torque at the base of the link and command control signals that adjust the motor references θ * m (t).

Estimators: (a)
The impact detector monitors the motor's position and force-torque at the base of the link to detect the instant at which the antenna impacts an object.(b) The contact point estimator determines the point of the antenna at which contact has been detected.
The combination of these controllers and the estimators across the three stages of the control process is described.
In the first stage (free motion control), the inner loop is activated along with the impact detector, which continuously monitors the data.A programmed motor trajectory allows the antenna to perform a sweep.Then, the antenna moves freely until it makes contact, at which point the impact detector activates and triggers the transition to the second stage.
During the second stage (post-impact), a new reference is set for the inner loop, causing the motors to increase their position relative to the angle at which the contact has been detected, thereby ensuring that the antenna continues to exert pressure on the object.Then, the system remains steady for a predetermined period of time, collecting force-torque data from the base of the link until the contact point estimator identifies the point at which the antenna is pushing against the object.
Following contact point estimation, the third stage (force control) commences.The force control strategy incorporates two nested control loops: the inner loop and the outer loop.The inner loop regulates the motor position, while the outer loop indirectly controls the exerted force by adjusting the torque at the base of the link.Once the distance from the contact point to the joint is estimated, the outer loop sets the desired torque at the base of the link as a reference.This reference is calculated by multiplying the estimated distance by the desired force.
Detailed descriptions of all subsystems involved in the control process are provided below.

Motor Control Inner Loop
The inner loop is designed to control the position angle of the actuators so that the dynamics between the motor position and its reference become an approximately linear time-invariant system.This control is insensitive to gravity and external forces acting on the antenna and remains active throughout the entire control process.This control system has been utilized in both free-and constrained-motion scenarios (e.g., [39]), demonstrating its robustness and effectiveness.The structure incorporates PID controllers with a lowpass filter term, ensuring excellent trajectory tracking, compensating for disturbances such as unmodeled friction components, and maintaining robustness against parameter uncertainties.This provides precise and rapid motor positioning responses.Additionally, it includes a compensator for the nonlinear friction of the motor (stiction) to avoid motor dead-zones and a compensator for the estimated coupling torque caused by the force exerted by the link on the contacted object.
An algebraic design methodology allows for the arbitrary placement of the four poles and two zeros of the closed-loop system.By aligning the zeros and poles at the same position, denoted as p m , and assuming that the compensators perfectly cancel the nonlinearities, the inner closed-loop dynamics are defined as follows: This configuration allows for very rapid motor movements if the absolute values of the poles p m are set high, provided that actuator saturation is avoided.

Impact Detector
Determining the precise moment of contact is crucial for initiating the mechanisms that estimate the contact point on the beam.In robotics, several mechanisms have been proposed that detect a collision by monitoring measured variables that exceed a threshold: [40] for rigid-link robots and [41] for flexible-link robots.Contact instants have also been estimated in artificial antennae designed to mimic insect behavior [42].These experiments employed a two-axis acceleration sensor positioned at the antenna's tip to measure link vibrations and gather information about contact.Information regarding contact was derived from analyzing vibration frequencies.In the case of our flexible antenna, we utilized a mechanism that predicts real-time coupling torque during the antenna's free movement.This mechanism involves comparing the predicted torque with sensor measurements in real time.The equation used to estimate coupling torque is derived from the dynamics of the flexible antenna link.A brief outline of this estimator is given below.

Consider the measured coupling torque
⊺ provided by the F-T sensor at the base of the link.Denote the effect of gravity on the beam provided in the F-T sensor frame as − → Γ g (t) = (0, 0.5 being the mass of the antenna and θ m 2 (t) the elevation motor angle.Define − → Γ e (t) as a real-time estimation of the coupling torque during free-movement mode, assuming no gravity, obtained from Equations ( 31) and (32).The residual error between the measured and the estimated coupling torques can be defined as: Then, contact is produced at the instant t i at which the absolute value of the time derivative of the magnitude of the residual error exceeds the threshold: where the threshold r max Γ is determined experimentally.

Contact Point Estimation
We propose determining the contact point of the antenna, where it makes contact with the surface of an object, using the algorithm outlined in [39].This algorithm combines two estimators.The first one relies on the relation between the lowest natural frequency of the oscillations experienced by the antenna after an impact, ω 1 , and the contact position, l c , as described in Section 4 and represented in Figure 4.This relationship is usually tabulated, allowing for a quick estimation of the contact point from the frequency value.Although this method gives a very precise estimation, it can sometimes yield two possible solutions.To resolve this ambiguity, a second estimator is employed.This one estimates the contact point using the static force and torque measurements of the sensor and the relation between these magnitudes.Since torque is the product of force and distance, it straightforwardly determines the application point.While this method may be less precise than the first, it effectively distinguishes between the two potential solutions provided by the initial estimator.
The contact point estimation process begins when the impact detector triggers the transition from the first stage (free motion of the antenna) to the second stage (post-impact).The antenna pushes against the object and remains steady for a determined period of time ∆t, during which the F-T sensor registers the oscillations of the antenna.Subsequently, Fast Fourier Transform (FFT) is performed on the data to determine the first vibration frequency ω 1 , and the contact point l c is obtained from the tabulated data, as depicted in Figure 4.In cases where two potential contact points are identified, the second estimator computes the contact point, determining which of the initial estimates is correct.
The precision of obtaining the frequency ω 1 by performing FFT on the registered data is inversely proportional to the length of the data vector.Longer data vectors provide greater precision, but also increase algorithm execution time.Therefore, a balance must be struck in defining the data registration time, ∆t.

Force Control Outer Loop
Force control is indirectly achieved by controlling the torque at the base of the antenna.If the force we aim to exert is F * (t) at contact point l c , then a moment at the base of the antenna Γ * (t) = F * (t) • l c should be exerted.The feedback measure from the F-T sensor Γ(t) and the rest of the measures are modified to correct the gravity effect of the antenna system, as is performed for the impact detector in Section 6.2.Thus, the feedback signal is obtained from ⊺ being the measured coupling torques and forces, respectively, provided by the F-T sensor at the base of the link, and − → Γ g (t), − → F g (t) the torque and force effects of gravity on the beam provided in the F-T sensor frame, where the product ρ • L is the mass of the antenna and θ m 2 (t) the elevation motor angle.
Figure 6 presents the outer loop of the system, which is composed of the following: 1.The transfer function G * (s, λ c , θ e ) from Equation (77) describing the dynamics of the antenna in contact with an object.2. The motor control inner loop G M (s) whose dynamics are described by Equation (82).3. A controller C(s) whose robustness would be justified if it were of the PI type: which verifies the robustness condition Figure 6.Force control outer loop scheme.
This control system robustly stabilizes the dynamics of the robot under contact, maintaining stability for any pair of values where 0 < λ c ≤ 1 and −90 • ≤ θ c ≤ 90 • , as well as for any uncertainties in the mechanical parameters of the antenna.The next subsection will prove the above robustness conditions.

Stability Robustness Condition of C(s)
The robust stability of the PI force controller (87) is designed using the Routh-Hurwitz criterion, e.g., [43].
The closed-loop characteristic equation is We hereafter omit the arguments of K a , β, and α for the sake of simplicity.Then, the characteristic polynomial is In order to assess the closed-loop stability, we must first check that all the coefficients of this polynomial are positive.This is easily verified, since α, β, K a , K c , a c , and ε are positive.Next, we calculate the Routh table, giving the following terms in the first column: It is easy to see that 2 − a c • ε > 0 with α > β are sufficient conditions to make this term positive.

•
Term where And, again, it is easy to see that 2 − a c • ε > 0 is a sufficient condition to make this term positive.
Therefore, considering that the property α > β seen in the previous section is verified, we have proven that a PI controller provides force control with robust stability if the following condition 0 is verified.

Design Methodology of C(s)
In this particular application, an algebraic methodology is followed to adjust the parameters K c , a c of the controller C(s) (87).The characteristic equation of the system is Calculating Equations ( 93) with (87), ( 82) and (78), while omitting λ c , θ e for clarity, yields: The two parameters K c and a c of the controller C(s) need to be adjusted.Thus, if a double pole of the system is selected and placed at p F , the following relations need to be accomplished: From Equation (95), the parameter K c based on a c is obtained: Finally, by calculating Equation (96) and including (97), the parameter a c is obtained:

.3. Justification on Tuning the Outer Loop Considering the above Robust Stability Condition
Expressions (97) and (98) enable real-time tuning of the PI controller (87) once λ c and θ e have been estimated.We recall that θ e has minimal influence on the parameters of transfer functions G a (s, λ c , θ e ).However, if necessary, it can be determined by combining measurements from the motor encoders and an inclinometer mounted on the base of the haptic system.The tuning process aims to achieve the same closed-loop poles indepen-dently of the contact point λ c .Nevertheless, this ideal situation cannot be fully achieved in practice due to the following factors: 1. Variations in the estimation of the frequency ω 1 from measured data can lead to incorrect contact point estimations.As previously mentioned, the precision of the FFT depends on the length of the data vector, defined by ∆t.Since this period cannot be set too high, it inevitably introduces imprecision in contact point estimation.2. Modelling errors cause variations in the curve relating the contact point and first vibration frequency of the antenna, represented in Figure 4. 3. The transfer function G a (s, λ c , θ e ) (77) obtained from the model is a simplification of the full system, truncated to the first mode of vibration.This simplification introduces modeling errors in the parameters of this transfer function, leading to non-optimal calculations of controller parameters.
These issues can result in unstable outer loop control systems if a robustness condition is not imposed in the design of C(s).An unstable outer loop can cause undesirable and dangerous behavior, potentially exerting excessive force and risking the integrity of the robot.We have proven that condition (92) guarantees closed-loop stability in all these cases.Moreover, it ensures limited deterioration of the transient response in cases of mismatch.

Robot Parameters and Experimental Results
Experiments are conducted to test various contact points and programmed reference forces.Both degrees of freedom of the robot (azimuthal and elevation movements) are evaluated.Specifically, a set of seven different contact points ranging from λ c = 0.3 to λ c = 0.9 and three levels of reference forces |F * | = (0.05, 0.10, 0.15) N are tested between one and five times for each degree of freedom.An image of the experimental setup is shown in Figure 7, where the system is performing an azimuthal (horizontal) movement with the sensing antenna making contact with a steel cylinder at λ c = 0.9 of the antenna.
In this section, we first detail the main parameters of the system and the control process.Then, we present the experimental results obtained from the different algorithms at each stage of the control process.In this photo, the system is performing an azimuthal (horizontal) movement with the sensing antenna hitting the steel cylinder at λ c = 0.9 of the antenna.

Parameters of the System
Table 1 shows the parameters of the two motors of the system, where Motor 1 and Motor 2 refer to the azimuthal and elevation motors, respectively.
Table 2 details the characteristics of the antenna.Note that the link flexural rigidity EI, as defined previously, is a product of Young's Modulus E and the area moment of inertia I. Finally, Table 3 shows the most important parameters of each control system: • First, in the motor control inner loop, the closed-loop system's poles p m are placed at the same value for both the azimuthal and elevation motors to achieve homogeneous behavior of the system in both degrees of freedom.• Second, the impact detector threshold r max Γ , which has power units (N • m/s), is determined experimentally based on the maximum value of (84) obtained in the free-motion experiments, with an added security margin.• Third, in the contact point estimator, a time of ∆t = 0.7 s is chosen as it provides sufficient FFT precision while allowing the algorithm to execute quickly enough.In this case, the relation between ω 1 and l c , described in Section 4 and represented in Figure 4, is tabulated to allow for quick estimation of the contact point from the frequency value.• And fourth, in the force control outer loop, the closed-loop system's poles p F are placed to achieve the fastest outer loop response possible while satisfying the robustness condition (92).Furthermore, the different values of the parameters K c (λ c , θ e ) (97) and a c (λ c , θ e ) (98) of the force controller C(s) (87) are tabulated to facilitate quick tuning of the outer loop during experimentation.The results of the experiment depicted in the photo in Figure 7, where the antenna performs azimuthal displacement, contacts the cylinder at λ c = 0.9, and pushes with a programmed force of F * = 0.15 N, are represented in Figure 8.The data illustrate the complete control process, from the first stage of free motion control, through post-impact data acquisition in the second stage, to the force control in the third stage.Hereafter, the graphical results presented below belong to this same experiment.

First Stage Results: Impact Detector
The time required for the impact detector to detect contact is measured using a special setup involving the object that the antenna impacts.This setup consist of a thin copper wire attached very close to the surface of the steel cylinder, but not touching it.The cylinder is wired to the digital input of the DAQ system and is set to zero volts.The wire is connected to an output port of the DAQ supplying 5 volts.When the antenna hits the cylinder, it also pushes the wire towards the cylinder surface, causing an electrical connection between the wire and the cylinder.This results in a voltage change in the digital input of the DAQ system connected to the cylinder, registering the exact instant t i A at which the antenna contacts the cylinder.This setup is hereafter referred to as the analog impact detector, and a detailed image of it is shown in Figure 9. Figure 10 shows the performance of the system during the first stage, where the motors control the movement of the antenna, seeking the space until it makes contact with the cylinder.The figure includes plots of the motor reference versus encoder signal, measured torque versus torque simulated by the detector, and measured force applied to the cylinder.It also represents the moment at which the analog impact detector (t i A ) detects contact.Table 4 summarizes the mean results and the standard deviation of all experiments regarding the delay in estimating the contact instant.The time required for the impact detector to detect contact is calculated as ∆t i = t i − t i A , where t i > t i A .Alongside this table, a histogram of ∆t i for all experiments is shown in Figure 11.The histogram illustrates that the most frequent time estimation delay falls between 0 and 2 milliseconds.The experimental setup is positioned in each experiment for the impact to occur at a programmed position of the antenna.Specifically, a set of seven different contact points from λ c = 0.3 to λ c = 0.9 are measured and marked on the antenna with a white point (see Figure 7).Figure 12 shows the collected data during this second stage, where the antenna remains steady, pushing the cylinder for a determined period of time ∆t.This parameter determines the precision with which the FFT determines the frequency.The precision in the FFT frequency is calculated as the maximum frequency read f max , which is half of the frequency of system, divided by the length of the registered data L data , related to ∆t such that Taking into account the relation between ω 1 and l c described in Section 4 and represented in Figure 4, the maximum ∆ f that can be selected varies between 1.4 and 1.5 Hz.This data are obtained considering less than a 2% error in the length of the antenna when estimating the contact point, which corresponds to approximately 10 mm.Thus, the data acquisition time selected for the second stage is ∆t = 0.7 s.Table 5 summarizes the mean values of the estimated contact point and its errors for all the experiments.It can be observed that the mean absolute errors do not exceed the limit of 10 mm, which is approximately 2%.Alongside the table, Figure 13 shows the estimated contact points l c for all experiments in comparison with the real contact point reference l * c = λ c • L. It can be seen that the estimator provides accurate results.14 illustrates the performance of the control system during this third stage, where the antenna applies a specified force to the cylinder.During this phase, both the inner and outer loops of the control system operate concurrently.The programmed force is the reference of the outer loop, and the control signal that it generates results in the input reference of the inner loop.Figure 14 demonstrates that the system operates effectively with zero steady-state errors.
As explained earlier, the outer loop controls the torque at the base of the antenna with the reference Γ * (t) = F * (t) • l c , where F * (t) is the desired force and l c is the estimated contact point from the previous stage.A PI controller is used, which achieves zero steadystate error in torque for every experiment conducted.However, variations in the estimation of the contact point l c lead to two issues: (1) tuning C(s) with non-optimal parameters, and (2) setting an incorrect reference torque Γ * (t).The first issue affects the transient response of the resultant system, as it does not operate as quickly as theoretically predicted.Ideally, the settling time t s of the outer loop response, obtained from simulations, is t s = 0.135 s.This ideal result can be compared with the data in Table 6, which presents the mean settling time measured in each experiment.Additionally, Figure 15 shows a histogram of the settling times t s obtained for all experiments.Finally, the second issue affects the steady-state response of the system.Since the estimation l c of the contact point may introduce errors compared with the real contact point l * c (see Figure 13), setting the reference torque Γ * (t) causes the control to push the cylinder with a force of F(t) = F * (t) • (l c /l * c ), which does not exactly correspond to the desired applied force |F * | = (0.05, 0.10, 0.15) N. The percentage error between the desired force F * and the real applied force F is calculated as: The mean percentage error of force is calculated considering the errors introduced by the estimation of the contact point in all the experiments conducted, and it is illustrated in Figure 16.These results are consistent with the errors observed experimentally (e.g., in Figure 14).

Discussion
This paper developed and tested a precise force control mechanism for a haptic device comprising a flexible link that rotates around one of its ends, resembling the antennae found in many insects.The link, with a distributed mass, executes azimuthal and elevation movements influenced by gravity.Contact with an object can occur at any intermediate point along the link.It is crucial in this context to regulate the force exerted by the antenna on the object to facilitate tasks such as object identification or moving an object.Previous works have addressed force control only when the contact is at the tip.Our work makes several significant contributions to the state of the art because, for the first time, (1) precise force control at intermediate points of a link is achieved; (2) a condition to design robustly stable controllers is obtained, i.e., controllers that maintain acceptable performance independently of the features of the controlled dynamics, that highly change with the contact point at the link; (3) we prove that simple PI controllers verifying this condition achieve such robustness stability; and (4) this control system yields satisfactory experimental results.Moreover, a lumped-mass model (with more than a lumped mass) of a flexible link in the context of contact with an object at an intermediate point is developed for the first time.This model is general because it is developed for a normalized beam.
Next, we specify the roles played by the dynamic models developed in this work.In the first scenario, in which the antenna moves freely, vibrating without suffering any contact, the obtained free model is used to predict (by simulating this model) the coupling torque.This prediction serves to compute the residue used by the impact detector for estimating the impact instant.In the second scenario, the antenna presses the object in a motor control fashion, without employing a flexible-link model.In the third scenario, the parameters of the PI controller are tuned using the family of approximate models obtained for the case of contact (these contact models are different in terms of the function of λ c , and the forms of their dependence on λ c vary in terms of the function of the interval between masses that is being considered or whether λ c corresponds to the position of one of the masses of the lumped dynamics model).Moreover, we mention that the obtained contact models played decisive roles in obtaining the family of robust controllers: (1) these models yielded truncated models with two imaginary poles and two imaginary zeros that were used in the closed-loop stability assessment, and (2) they allowed us to establish a property whereby the zeros are closer to the origin of the S-complex plane than the poles, a property that was crucial for obtaining the robustness condition.
We highlight that we have designed a broader control system in which the PI force controller is embedded.It includes also an impact detector and a real-time estimator of the contact point.The experiments conducted using this whole system demonstrate the effectiveness of this methodology, ensuring the stability of the system and achieving minimal force error at the contact point.Figure 14 shows a mean error in the steady state of nearly null as consequence of using a PI controller.However, since the values of the exerted forces are low, noticeable noise can be observed in the figure because these values are not far away from the accuracy level of the force-torque sensor.
Next, we mention some limitations of the system.The first one is the precision of the controlled force, which depends on the strain gauge offset of the F-T sensor and, as previously mentioned, on the inaccuracy produced in the estimation of the contact position, which introduces a small error in the calculated torque reference.Another limitation is the assumption of small deflection.If this assumption were violated, the model obtained in Section 3 would be incorrect and the dynamics would become nonlinear.Finally, a third limitation is the assumption of the constant cross-section of the antenna.Other behaviors can be obtained using conical antennae.In this case, Section 3 should be developed assuming a decreasing cross section radius.
Finally, we mention that potential applications of the proposed force control exceed the haptic antennae case.This can be applied in other robotic scenarios like the following: (1) in biomimetics, where it can be used to design robotic birds that grasp objects with their beaks; (2) in industrial robots, where it can be used to design hands with flexible fingers that grasp objects with a programmed force, with contact at intermediate points of the fingers; and (3) in robot-assisted surgery, where a required force has to be exerted when the robot contacts an organ.quency, shown in Figure A1a, and this frequency tends to infinity as the contact approaches the tip of the beam, i.e., the position of the mass.In the case of the model with n = 2 masses, the minimization of (66) is performed by varying µ 1 and λ µ 1 in the interval (0, 1) with a step of 0.01.This gives a minimum mean square error of MSE = 0.390 for the model whose parameters are µ 1 = 0.63, µ 2 = 0.37, λ µ 1 = 0.39, and λ µ 2 = 1.However, a local minimum with MSE = 0.410 is highlighted for the model with parameters µ 1 = 0.63, µ 2 = 0.37, λ µ 1 = 0.72, and λ µ 2 = 1 whose MSE is close to the global minimum.These two models have two vibrational frequencies, which are shown in Figure A1a,b.In this case, it is the second frequency that tends to infinity when the contact coincides with one of the masses.
Figure A1 shows the three frequencies of the model and, as in the previous cases, the third frequency tends to infinity when contact occurs at the position of one of the masses.

Figure 1 .
Figure 1.Haptic sensors of the antenna and whiskers types in nature.

Figure 3 .
Figure 3. Scheme of the flexible beam (a) prior to contact and (b) after contact.

Figure 7 .
Figure 7. Experimental setup.In this photo, the system is performing an azimuthal (horizontal) movement with the sensing antenna hitting the steel cylinder at λ c = 0.9 of the antenna.

Figure 8 .
Figure 8. Complete experimental results.Motor angular position (inner loop), measured torque and force (outer loop) along the three stages of the control process.Case: azimuthal displacement, contact point λ c = 0.9, programmed force of F * = 0.15 N.

Figure 16 .
Figure 16.Summary of stage 3 results (II): mean percentage force error obtained for azimuthal and elevation experiments.

Figure A1 .
Figure A1.Vibration frequencies as a function of the contact point, where (a) shows the first frequency, (b) the second, and (c) the third.Here () corresponds to the model of[25], ( ) to the model with n = 1 masses, () and ( ) to n = 2 masses for λ µ 1 = 0.27 and λ µ 1 = 0.72 respectively, and ( ) to the model with n = 3 masses.

Table 1 .
Parameters of the motors.

Table 2 .
Parameters of the antenna.

Table 3 .
Parameters of the control system.

Table 4 .
Summary of stage 1 results: mean and standard deviation σ of all experiments regarding the delay in estimating the contact instant.Histogram of the delay in estimating the contact instant ∆t i for all experiments.7.2.2.Second Stage Results: Contact Point Estimation

Table 5 .
Summary of stage 2 results: mean values of the estimated contact points (in millimeters) and their errors (in millimeters and % with respect to L) for all the experiments.

Table 6 .
Summary of stage 3 results (I): mean of the settling time responses t s of the outer loop in seconds.Histogram of the settling time responses t s of the outer loop.

Table A1 .
Root mean square error of fitted functions.